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Abstract 

The coefficient ca required for 0(a) improvement of the axial current in lattice QCD with 
A r f = 3 flavors of Wilson fermions and the tree-level Symanzik-improved gauge action is 
determined non-perturbatively. The standard improvement condition using Schrodinger 
functional boundary conditions is employed at constant physics for a range of couplings 
relevant for simulations at lattice spacings of ~ 0.09 frn and below. We define the improve¬ 
ment condition projected onto the zero topological charge sector of the theory, in order 
to avoid the problem of possibly insufficient tunneling between topological sectors in our 
simulations at the smallest bare coupling. An interpolation formula for ca(#o) is provided 
together with our final results. 
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1 Introduction 


Wilson fermions [1] are an attractive and popular fermion discretization for lattice QCD 
simulations. However, while the addition of the Wilson term lifts unwanted fermion ‘dou¬ 
bler’ modes, it also leads to 0(a) cutoff effects. At small a, these cutoff effects are described 
by the (continuum) Symanzik effective theory [2,3], and the procedure of their systematic 
removal is called ‘0(a) improvement’. 

In order to eliminate these 0(a) cutoff effects, a single dimension-five ‘clover’ term 
must be added to the Lagrangian [4] (with coefficient c sw ), and additional counterterms 
have to be added to local operators. To this end, the (unrenormalized) improved axial 
current is given by 


( A i )£ = A«(x) + ac A ^(d fl 

+ d;)P a (x), 

(1.1) 

A*(x)=$(x)T a 7m 75^(x), 

P a = ^(x)T a 75 ^(x), 

(1.2) 


where T a is an SU(IVf) generator acting in flavor space and d^, d* are the lattice forward 
and backward derivatives, respectively. Therefore, in order to improve matrix elements of 
the axial current, the coefficient ca must be specified as well as c sw . The non-perturbative 
determination of ca for a range of bare couplings is the subject of this work. 

Matrix elements of the axial current are of particular importance, since they enter the 
computation of pseudoscalar meson decay constants and quark masses. These quantities 
are not only of great phenomenological interest in their own right, but it has also been 
demonstrated [5] that the kaon decay constant /k provides a precise determination of the 
lattice scale in physical units, thus affecting all dimensionful quantities. For instance, in the 
two-flavor theory and at lattice spacings of ~ 0.1 fm, the CA-ternr discussed above typically 
contributes ~ 10 — 15% to pseudoscalar meson decay constants and quark masses [6,7]. 

In this work we consider lattice QCD with iVf = 3 mass-degenerate flavors of Wil¬ 
son fermions and the tree-level Symanzik-improved gauge action [8]. This gauge action 
is demonstrably preferable in the pure gauge theory, where its cutoff effects are smaller 
than other standard actions [9]. In preparation for dynamical simulations in this setup 
(see [10] for a first report), the parameter c sw multiplying the dimension-five ‘clover’ term 
in the action has been tuned non-perturbatively [11]. Like c aw , the improvement coefficient 
ca has been determined in perturbation theory at 1-loop in Ref. [12], However, previous 
non-perturbative determinations of ca [6, 7] deviate roughly 300 — 400% from 1-loop per¬ 
turbation theory at the largest bare couplings, so that a non-perturbative determination 
is required. 

To determine ca non-perturbatively, we employ the (by now) standard improvement 
condition for dynamical fermions [6, 7], which imposes the PCAC relation at constant 
physics to ensure a removal of O(a) effects in on-shell quantities and, at the same time, 
a smooth behavior of vanishing 0(a 2 ) effects as the bare coupling is varied. However, for 
technical reasons related to the topology freezing of our simulations at the finest lattice 
spacing, we project the necessary correlation functions onto the trivial topological sector. 
The main result of this work is the interpolation formula for ca(5o) gi ven hi Eq. (4.1), with 
the coefficients of Eq. (4.2), which is valid for lattice spacings of a ~ 0.09 fm and below. A 
statistical error of 4% should be assigned to this formula near the largest simulated bare 
couplings, whereas 7% is more appropriate at the finest lattice spacing. 
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We detail the improvement condition as well as the projection onto the trivial topo¬ 
logical sector in Sect. 2, and our simulation setup is discussed in Sect. 3. Numerical results 
together with the final interpolation formula are presented in Sect. 4, and conclusions are 
drawn in Sect. 5. 

2 Improvement condition 

We determine c\ via a variant of the improvement condition originally introduced in 
quenched QCD in [13] and applied in its present form to the theory with dynamical fermions 
first in [6] for the two-flavor case. We briefly review this condition here, adopting the no¬ 
tation of Ref. [6]. 

Improvement conditions are typically based on imposing the PCAC relation, which in 
its continuum form is an operator identity, at finite lattice spacing. A consequence of this 
relation is that the PCAC quark mass, defined as 

” P0AC = 2{a\P°( x )\l3) ' 

is independent of the space-time position x as well as the external states a and (3. It is 
this property of mpcACo required to hold on the lattice up to 0(a 2 ) cutoff effects, which 
we exploit to determine c\. 

Employing the above relation with the improved axial current of Eq. (1.1), we can 
decompose 


( 2 . 1 ) 


mpcAC = m(x; a, f3) = r(x; a , /3) + qca ■ s(x; a, /3), (2.2) 


r(x;a,(3) 


{<x\h( d n + d »)(A(x))%\p) 
2 (a \P(x) a \ (3) 


s(x; a, / 3) 


{g\d^{P(x)) a \p) 
2 (a \P(x) a \/3) 


(2.3) 


As mentioned above, if the continuum form of the PCAC relation holds, m(x; a, f3) = m 
is independent of a, (3 and x. If we demand this property and fix x, the quark mass m 
can be eliminated by using two different pairs of external states a, j3 and 7, 5 to obtain our 
definition of ca- 


1 r(x;a,P) — r(x;'y,6) 
a s(x-,a,f3) — s(x;'y,5) 


(2.4) 


In practice we employ a finite physical system size (L ~ 1.2 fm) with Schrodinger 
functional boundary conditions [14,15] in time and a periodic torus in space. Furthermore 
(using the standard notation), we employ a vanishing background gauge held (</>* = </>( = 0) 
and periodic spatial boundary conditions for the fermion fields (#* = 0). With the tree-level 
Symanzik-improved gauge action, there are several possibilities for implementing boundary 
0(a) improvement in the Schrodinger functional. As in Ref. [11], we resort to ‘choice B’ of 
Ref. [12], which possesses the desirable property that the classical minimum of the action 
can be expressed analytically. Our boundary 0(a) improvement is implemented at tree- 
level according to this choice. An additional boundary 0(a) improvement possibility for 
this action in the Schrodinger functional has been proposed recently in Ref. [16]. 
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Correlation functions in the Schrodinger functional setup may involve boundary quark 
fields, which are adopted here. We construct r and s defined above from the following 
correlation functions: 

/a(w; U ) = -^X>S(x)0“M), /p(x 0 ^) = -^X;<P“(x)O>)), (2.5) 

X X 


O a (uj) =a 6 ^2 CO?) • T a 75 • u(x - y) • C (y), (2.6) 

S,y 

where the pseudoscalar operator O a is composed of the boundary quark fields at xo = 0 
(C and C), while the ‘wavefunction’ c o(x) will be discussed shortly. Dehning the correspond¬ 
ing operator 0 /a (u/) at xq = T, we also employ the boundary-to-boundary correlation 
function 

(2.7) 

As outlined before, we aim at probing the PCAC relation with two different pairs of 
external states. To achieve this, we construct approximate wavefunctions of the ground 
and first excited state in the pseudoscalar channel, viz. 

3 3 

w 7r(°) ~ 5 ^i 0 ) w», w w( i) » ( 2 . 8 ) 

i— 1 2—1 

where 0J n (o) (uJ n ( i>) is constructed to maximize the overlap with the ground (first excited) 
state. These approximate wavefunctions are superpositions of trial wavefunctions with 
coefficients 77 ^, given by the eigenvectors of the matrix / ? y = f\ corresponding to 

the largest and next-to-largest eigenvalues for w ( 0 ) and w (i), respectively. As our spatial 
trial (hydrogen-like) wavefunctions at the boundaries we take 

u>i(r) = e - r / ao , tc 2 (r) = r-e- r / a °, Co 3 {r) = e~ r/ ^ ao) , (2.9) 


!Xi(x) = N t Wj (\x - nL I), 

nEZ 3 


( 2 . 10 ) 


where N t is a normalization factor chosen to ensure a 3 ^Wu;?(x) = 1 and ao = L/6. 
Finally, our operational definition of ca now reads 


ca{xq) 


1 _ r(xp; io,Lo n (i)) - r(x 0 ; i 0 , u wi0) ) 

a s(x 0 -,io,u) n (i)) — s(x 0 -,io,u) n ( 0 )) 


1 Ar(xo) 
a As(xq) ’ 


( 2 . 11 ) 


r(x;i 0 , u) 


l{do + d^)f A (x 0 -,oj) 
2/ P (x 0 ;w) 


s(x;z 0 ,w) 


dodo/pQro;^) 

2/p(x 0 ;w) 


( 2 . 12 ) 


which at the same time defines Ar(xo) and As(xo). In this way, the wavefunctions cj ( 0 ) 
(u}^(i )) determine the states j3 ( 5 ), cf. Eq. (2.4), while for both states a and 7 the plain 
Schrodinger functional boundary state i 0 with vacuum quantum numbers is inserted. The 
choice of xq will be discussed later. 
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Unlike previous studies, we employ T = 3L/2 lattices for the generation of our dynam¬ 
ical gauge field ensembles, with the intention of re-using them for a determination of the 
axial current renormalization constant Za according to the method of Ref. [17] (see [18] 
for a preliminary report). While in a quark mass independent improvement scheme as 
applied here one ideally would impose the improvement condition of Eq. (2.11) at zero 
quark mass, in practice it was found in Ref. [6] and is also confirmed by our data that ca 
is rather insensitive to fairly small deviations from this constraint. However, that is not 
the case for Za, so in general we endeavor to tune the quark mass to values close to zero. 

Like periodic temporal boundary conditions, Schrodinger functional boundary con¬ 
ditions are ‘closed’, and disconnected topological sectors emerge in the continuum limit. 
However, for small physical volumes, non-trivial topological sectors receive a small weight 
in the partition sum. Unfortunately, to keep O(a 2 ) effects under control, and with an 
eye toward using our ensembles for Za, our physical volume (L ~ 1.2 fm) is large enough 
that non-trivial topological sectors are not completely suppressed. Whereas these sectors 
are sampled sufficiently at coarser lattice spacings, ensembles at our finest lattice spacing 
(L/a = 24, f3 = 3.81) are effectively frozen in the sector with topological charge Q = 0. 

The issue of topology freezing in the Schrodinger functional has been investigated 
recently in the pure gauge theory in Ref. [19]. There it was suggested that quantities pro¬ 
jected to the zero topological sector have a smooth approach to the continuum limit. For 
the case at hand (ca as well as Za), quantities defined in this way differ from their conven¬ 
tional counterparts only by irrelevant cutoff effects. Formally, we perform this projection 
for all observables of interest. Adopting the notation of Ref. [19], we define 

(O )o = , (2.13) 

\°Q, o) 

where O is an arbitrary observable in the theory (such as the correlation functions in 
Eqs. (2.5) and (2.7)), and the topological charge Q is defined using the Wilson (gradient) 
flow [20, 21] at flow time t given by y/8t/L = c with c = 0.35. Since at finite lattice 
spacing the topological charge may take non-integer values, we define all configurations 
with \Q\ < 0.5 as the trivial topological sector. Thus (as in Ref. [19]) we make the 
replacement 5q : q —> 9{Q + 0.5)0(0.5 — Q). For ease of notation we shall henceforth take 
ca to mean the one projected onto the trivial sector in this manner, with the exception of 
Tab. 2, which directly compares results in all sectors (where available) with those restricted 
to Q = 0. Let us anticipate already here that these two kinds of analyses yield consistent 
results for ca as expected, because the Ward identities underlying the PCAC relation and 
thereby our improvement strategy hold in any topological sector. 

An alternative to address topology freezing in the Schrodinger functional has recently 
appeared [16], namely the use of ‘half-open’ boundary conditions. While not considered 
here, these boundary conditions may help the problem in future calculations provided an 
improvement condition which does not require boundary-to-boundary correlation functions 
is devised. 
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L 3 x T/a 4 

P 

Av 

#REP #MDU 

ID 

12 3 x 17 

3.3 

0.13652 

10 

10240 

Alkl 



0.13660 

10 

13048 

Alk2 

12 3 x 19 

3.3 

0.13652 

10 

10468 

A2kl 

16 3 x 23 

3.512 

0.13700 

2 

20480 

Blkl 



0.13703 

1 

8192 

Blk2 



0.13710 

3 

24560 

Blk3 

16 3 x 23 

3.47 

0.13700 

1 

8176 

B2kl 

20 3 x 29 

3.676 

0.13680 

1 

7848 

Clkl 



0.13700 

4 

15232 

Clk 2 



0.13719 

4 

15472 

Clk3 

24 3 x 35 

3.810 

0.13712 

7 

15448 

Dlkl 


Table 1: Summary of simulation parameters, number of replica and total number of molec¬ 
ular dynamics units of our gauge configuration ensembles labeled by ‘ID’. 

3 Simulation details 

Our simulations are performed using Schrodinger functional boundary conditions. We use 
the openQCD code 1 of Ref. [22], which was also used for the simulations to calculate c sw in 
Ref. [11]. 

The bare gauge couplings are chosen to approximately satisfy a constant physics condi¬ 
tion, fixing L ~ 1.2 fm. In this way it is ensured that any 0(a) ambiguities in c\ disappear 
smoothly toward the continuum limit. For a thorough an more general discussion of the 
idea and virtues of imposing improvement (and renormalization) conditions at constant 
physics, see, e.g., Refs. [23,24], Similarly to earlier work [6,7], we fix the physical vol¬ 
ume by beginning with a particular pair of and L/a (/3 = 6 /</q = 3.3 at L/a = 12 
in the present case) and choose the bare couplings for subsequent smaller lattice spacings 
according to the universal 2-loop /3-function 2 . The range of lattice spacings covered in 
this way extends from a ~ 0.09 fm to a ~ 0.045 fm. As already mentioned, at each bare 
coupling we tune the bare quark mass so that the PCAC mass is kept approximately con¬ 
stant and close to zero 3 . At several lattice spacings we confirm that our determination of 
ca is insensitive to variations of the (small) quark mass. Information about our ensembles, 
consisting of several replica per parameter set in most cases, can be found in Tab. 1. Due 
to practical reasons discussed in the openQCD documentation 4 , our lattices have temporal 
extents T = 3L/2 — a. Since this offset itself scales with a, one expects its influence on the 

1 http://luscher.web.cern.ch/luscher/openQCD/ 

2 Note that the non-universal 3-loop term of the /1-function is not known for the tree-level Symanzik- 
improved gauge action. 

3 Based on the experience from the two-flavor theory, for which the multiplicative quark mass renormal¬ 
ization factor only varies slowly with a, we can safely neglect it for the tuning purposes here, too. 

4 For this work we employ openQCD version 1.2. This issue has been corrected in the latest version (1.4). 
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determination of ca to be of O(a 2 ) and thus to be small; still, we assess it on our coarsest 
lattice where it is largest. There we simulate 12 3 x 17 as well as 12 3 x 19 ensembles. 

Although it was found previously that small deviations from the constant physics 
condition have little effect on ca [6 ], we estimate this deviation by measuring the scale- 
dependent renormalized coupling </q F , defined in Ref. [25]. Results for this coupling are 
shown in Tab. 2. We also test the dependence of ca on L in physical units (and thereby 
on violations of the constant physics condition) directly by simulating an additional bare 
coupling at L/a= 16. 

We now briefly summarize the simulation algorithm used for the generation of these 
gauge held ensembles. While two of the (mass-degenerate) pseudo-fermion fields can be 
simulated in the usual way, the RHMC algorithm [26] is employed for the third. Even-odd 
preconditioning is used for all fermion determinants, whereas mass preconditioning [27] 
with two additional pseudo-fermion fields is used for the degenerate doublet. We use a 
hierarchical integration scheme [28], where the gauge force is integrated on the innermost 
level and the remaining fermion forces on the second level. For the L/a = 16, 24 and part 
of the L/a = 20 lattices, the lowest poles of the RHMC are integrated on a third level. 
For the two inner levels, a fourth-order OMF integrator [29] is used, while the third level 
(when present) uses a second-order OMF integrator. A single step is used for the inner 
integrators, and the number of steps for the outer level is tuned to achieve an acceptance 
rate of ~ 90%. For the coarsest L/a = 12 lattice ensembles Alk2 and A2kl, we adopt 
(type I) twisted mass reweighting [30]. The conjugate gradient solver is employed for most 
fermion forces, while the multi-shift variant is typically used for most of the RHMC poles. 
For the lightest mass-preconditioned held and RHMC poles on the L/a = 16,24 lattices, 
we employ the SAP-preconditioned GCR algorithm [31]. 

4 Results 

On most of our ensembles, we measure the correlation functions defined in Eqs. (2.5) 
and (2.7) on every fourth trajectory of length r = 2MDU so that the spacing between 
these measurements is 8MDU. Only on Alk2 and A2kl, we use a measurement separation 
of 2 t = 4 MDU. The total statistics for all the ensembles considered here are tabulated in 
Tab. 1. 

In addition to these correlation functions, we also measure ‘smoothed’ gauge held 
observables obtained from the Wilson (gradient) how [20,21], which possess a well-defined 
continuum limit. These smoothed observables are useful in several ways. The smoothed 
gauge helds provide a renormalized definition of the topological charge, which we use to 
monitor the topology freezing discussed in Sect. 2. Even at lattice spacings where topology 
freezing is not a problem, the smoothed topological charge and action typically possess the 
largest observed autocorrelation times. Furthermore, the aforementioned renormalized 
(and L-dependent) coupling g^p of Ref. [25] is defined using the Wilson how and may be 
used to monitor the deviation from the constant physics condition, as it is sensitive to the 
physical lattice size. Results for this coupling are given in Tab. 2, too. 

In order to monitor the autocorrelation times in our simulations, we examine these 
smoothed observables at a how time t given by \/8 1/L = c with c = 0.35. For all simulations 
we find that integrated autocorrelation times of these observables satisfy the bound r max < 
200 —250 MDU, except for our L/a = 24 simulations where the charge is frozen. The other 


7 


.2 

ft 

« 


cd 

ft 

aj 


CD 

bjO 

!h 

cd 

ft 

CD 

23 

CD 

'5o 

o 

o 

ft 

o 



2200 

2000 

1800 

1600 

1400 

1200 

1000 

3 

2 

1 

0 

-1 

-2 

-3 

-4 


MDU 


P = 3.512, k — 0.13703 



0 4000 8000 12000 16000 


MDU 


s 

o 


4^ 

CD 

cd 

CD 

CD 

2 

ft 

'ft 


CD 

bO 

X 

ft! 

CD 

23 

CD 

'5o 

o 

o 

ft 

o 


1800 

1600 

1400 

1200 

1000 

1 

0 

-1 

-2 


1600 
1500 
1400 
1300 
1200 
1100 
1000 

1 

0 

-1 

-2 

0 2000 4000 6000 8000 10000 0 1000 2000 3000 4000 

MDU MDU 




Figure 1: Histories of the smoothed Wilson plaquette action and the topological charge for 
single representative replica from each of the ensembles, which enter into the final analysis. 
Our inability to sufficiently sample all topological sectors at /3 = 3.81 is evident. 
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ID 

ampcAC 

a?npcAC,o 

fi'GF 

fi'GF.O CA 

CA,0 

Alkl 

-0.0010(7) 

-0.0022(8) 

18.12(21) 17.77(20) -0.0551(26) 

-0.0594(31) 

Alk2 

-0.0086(6) 

-0.0100(8) 

16.95(13) 

16.62(15) -0.0557(19) 

-0.0552(24) 

A2kl 

-0.0011(7) 

-0.0025(10) 17.84(20) 

17.35(20) -0.0569(25) 

-0.0547(30) 

Blkl 

0.0063(2) 

0.0062(3) 

16.49(13) 

16.44(14) -0.0365(11) 

-0.0348(15) 

Blk2 

0.0056(3) 

0.0050(4) 

16.85(20) 

16.57(23) -0.0381(16) 

-0.0334(29) 

Blk3 

0 .0022(2) 

0.0016(3) 

16.11(14) 

15.78(15) -0.0380(11) 

- 0.0399(17) 

B2kl 

0.0041(4) 

0.0036(5) 

18.03(23) 17.95(26) -0.0342(22) 

-0.0344(46) 

Clkl 

0.0138(1) 

0.0137(2) 

16.55(27) 

16.44(26) -0.0324(14) 

-0.0305(29) 

Clk2 

0.0066(2) 

0.0065(3) 

15.53(14) 

15.40(15) -0.0300(21) 

-0.0311(26) 

Clk3 

-0.0005(1) 

-0.0006(2) 

14.64(13) 

14.41(16) -0.0281(14) 

-0.0291(14) 

Dlkl 

n.q. 

-0.00269(8) 

n.q. 

13.90(11) n.q. 

-0.0212(15) 


Table 2: Summary of results for ca- The (unrenormalized) PCAC quark mass ampcAC is 
computed from the correlation functions projected to the approximate ground state, using 
the 1-loop result for ca(9o) f rom [12], while g^ F denotes the gradient (resp. Wilson) flow 
coupling mentioned in the text. Recall that quantities with the explicit subscript label ‘0’ 
here refer to results from the analysis restricted to the sector of vanishing topological charge, 
whereas in the text we loosely suppress the ‘O’. Numbers for ensemble Dlkl (L/a = 24) 
are not quoted (‘n.q.’) for the case of covering all charge sectors in the partition sum, 
because our simulations are not able to sufficiently sample all the sectors and a reliable 
error estimation is thus not possible. Results from ensembles in italics enter into the final 
interpolation formula for ca(<?o), Eq. (4.1). 


smoothed observables turn out to still possess autocorrelation times of comparable order 
of magnitude in these simulations. The existence of the charge freezing at our finest lattice 
spacing is qualitatively illustrated in Fig. 1. There it is seen that the autocorrelation time 
of the smoothed action remains under control, whereas the autocorrelation time of the 
topological charge increases significantly from the L/a = 20 to the L/a = 24 ensembles. 
We are practically unable to sufficiently sample all topological sectors at this finest lattice 
spacing, necessitating the restriction of our observables to the trivial sector. 

To provide further support for this projection, we have estimated the expectation 
value {5q q); it effectively corresponds to the fraction of statistics, to which the restriction 
to Q = 0 reduces the number of generated configurations. For the ensembles with the 
smallest quark mass at given L/a, {Alkl, Blk3, Clk3, Dlkl}, we find (Sq.q) = 0.37(2) 
(L/a = 12), 0.44(2) (L/a = 16), 0.65(7) (L/a = 20), 0.90(5) (L/a = 24), and thereby 
to show large cutoff effects. These kind of cutoff effects have been observed before in the 
large-volume, two-flavor theory in Ref. [32], owing to the substantial suppression of the 
topological charge compared to the quenched approximation. In addition, our (5q } o) at 
m PCAC ~ 0 does not seem to smoothly (in the sense of monotonically and linearly or 
quadratically in a) approach 1 as L/a increases, which is the theoretical expectation for 
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(<5<2,o) in the large-volume continuum limit and at zero quark mass [33]. Even though we 
are not in a large-volume situation with mpcAC = 0 exactly, we interpret the encountered 
behavior of {5q i o) as a sampling problem of the algorithm, on top of the cutoff effects. 
Therefore, we prefer to project our results to the Q = 0 sector, which in the end does not 
induce a noticeable difference in the final numbers for ca, cf. Tab. 2. 

The estimation of the statistical error on our measured quantities is based on a single- 
elimination jackknife procedure, after first ‘binning’ the data (concatenated from different 
replica) such that the size of each bin is > r max , as well as on applying a full autocorrelation 
analysis according to Refs. [34,35]. Both yield similar error estimates; our quoted final 
results in the Q = 0 sector stem from the latter (without including any long-tail contribu¬ 
tions to the autocorrelation functions, which are negligible for the quantities entering the 
CA-analysis). 
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Figure 2: Top: Effective masses computed from fp with wavefunctions uJ n (o) and u f (i) 
(after projecting fp onto the approximate ground and first excited states) for the Blk2 
ensemble. The dotted horizontal line indicates the cutoff scale, L X a~ l = 16. Bottom: 
As(xo) and Ar(xo) for the same ensemble. 

After measuring the correlation functions, we solve for the largest two eigenvectors 
of the matrix given by f\ (ca( ,uij). These normalized eigenvectors have a well-defined 



- 1 1 1 1 1 1 1 1 1 1 




1 

) 

B 



B 




e 

- L/ a = 1 6 B m | I , 

I 1 

1 

1 

- 

0 = 3.512 





_ /c=0.13703 



- 

3 

T-H 

II 

c 

□ 

_1_L 


- 

. o n = 0: £J„ (0) 


- 

: e e e e e - - ® @ @ ‘ 


'_i_1_i_1_i_1_i_1_i_1_ 

_ 

' 


10 













continuum limit along our line of constant physics in parameter space, as long as the 
wavefunctions depend on physical scales only. In fact, as we do not observe any sig¬ 
nificant lattice spacing dependence for them, we fix these eigenvectors for once to the 
values calculated on the Blk2 ensemble (L/a = 16, (3 = 3.512, k = 0.13703) and 
regard them as part of our choice of improvement condition. For our setup we find: 
t/C°) = (0.5317(3), 0.5977(1), 0.6000(2)) and t/W = (0.843(5),-0.31(6),-0.44(6)), which 
are similar to those of Refs. [6, 7]. 

To get an idea of the sensitivity of our method to ca, we examine the effective masses 
of the correlation function fp, after taking the inner product with the eigenvectors for 
the (approximate) ground and first excited states. These are shown in Fig. 2 for the 
L/a = 16 ensemble Blk2 of the previous paragraph. The distinctly seen signals display 
that indeed the eigenvectors effectively maximize the overlap with the ground and first 
excited states, since these states are clearly separated up to xq ~ 11a. As already noted in 
Ref. [6], the energy of the first excited state is somewhat near the cutoff a -1 though, and 
this may influence the way in which residual cutoff effects are modified. Hence, residual 
0(a 2 ) effects may grow rapidly in smaller volumes [17], which justifies our choice of an 
intermediate volume with L ~ 1.2 fm to impose the improvement condition at constant 
physics. Also shown in Fig. 2 are Ar(.xo) and As(xo) (the latter being proportional to 
case °f exact ground and excited state projections) from Eq. (2.11) for the 
same ensemble, further demonstrating our good sensitivity to ca- 



X 0 /L 

Figure 3: The function ca(xq) f° r all the ensembles used in the final analysis. The dotted 
vertical line at xq/L = 0.5 indicates the space-time point used for our definition of ca- 

Finally, the function ca(xo) is shown in Fig. 3 for the ensembles used in the analysis. 
It indicates that ca is rather independent of xq for points sufficiently distant from the 
boundary. Moreover, only a rather little variation of ca(xo) is visible for x$ > 5a; this 
reveals that high-energy states, which could induce large O(a) ambiguities in the improve¬ 
ment condition, are reasonably suppressed in this region. As a compromise between cutoff 
effects and statistical errors, we take xq = L/2 as our final definition for ca, which besides 
is well within the regime where states with a distinct energy gap dominate the projected 
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correlators. These data are plotted in Fig. 5, together with an interpolation. 

Before discussing these final results for ca at each bare coupling, we assess our system¬ 
atic errors. In order to estimate the effect of our finite (but small) quark masses and thus 
of small violations of the constant quark mass condition on ca, we study several different 
sets of ensembles, which are identical apart from |Ltopcac| < 0.3 such as {Alkl, Alk2}, 
{Blkl, Blk2, Blk3}, and {Clkl, Clk2, Clk3}. Within each of these sets, the variation 
of the value of ca does not exceed more than about 1.5 standard deviations. To quantify 
the cutoff effect, which results from T = 3L/2 — a, we compare {Alkl, Alk2} with A2kl, 
where the temporal extent is T = 3L/2 + a. We see that this results in a difference which 
is significant at less than la at our coarsest lattice spacing where it is largest. Lastly, we 
assess the error due to the deviation from our constant physics condition. From Tab. 2 one 
infers that the ensembles {Alkl, Blk3, Clk3, Dlkl}, which enter into the final analysis, 
have < 7 q F 0 (L) between ~ 14—18, resulting in a ~ 20% variation assuming no cutoff effects. 
We can explicitly test the sensitivity of ca to this variation by means of the B2kl ensemble, 
which differs with respect to the B1 ensembles in a by ~ 6% and in <?q F 0 (L) by roughly the 
same amount as this overall g^ F 0 (L)-variation. The value of ca determined on the B2kl 
ensemble lies well within 1.5c of the interpolating formula in Fig. 5, so we are confident 
that this variation of L does not result in a significant shift of ca- Note that, after all, any 
imperfection in the constant physics condition of this level and the systematic uncertainty 
in ca induced by it will only introduce O(a 2 ) effects in quantities involving the improved 
axial current, which are negligible compared to other sources of errors and vanish in the 
continuum limit by definition. 



Figure 4: ca versus the lattice spacing measured in units of L for our ensembles {Alkl, 
Blk3, Clk3, Dlkl}. A linear fit to the data within the region where we simulate is also 
shown. 

We now move to the final results. In order to guide our choice for an interpolation 
formula, we observe that (for our choice of the improvement condition) ca is almost linear 
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Figure 5: Final results for ca together with the interpolation formula, Eqs. (4.1) and (4.2). 
The point from the B2kl ensemble is not included in the fit. Note that our non-perturbative 
results fall appreciably apart from 1-loop perturbation theory. 


in a over the range of bare couplings which we simulate. This approximate linearity is 
depicted in Fig. 4. Notice that the linear behavior can not extend all the way to a = 0, 
since this would be incompatible — owing to the non-polynomial relation between a and 
9o — with a polynomial dependence of c\ on in the perturbative regime. In our case, a 
naive linear fit to these data within the simulated region does not even extrapolate to zero, 
which is the value predicted by perturbation theory at tree-level: ca((7o = 0) = 0 [36,37]. 
On the contrary, we interpret the behavior of the results according to ca = constant + 
slope x a within the region of our data such that the (non-vanishing) constant term removes 
the targeted 0(a) effects in the non-perturbative regime, while the non-constant piece, 
describing the non-trivial dependence of c\ on g$, only affects 0(a 2 ) contributions in 
physical quantities. Regardless of their actual size, these intrinsic 0(a) ambiguities suggest 
to employ in the computation of physical quantities the -dependence of ca induced by 
the constant physics condition rather than, e.g., the constant piece alone, because one then 
expects the remaining leading 0(a 2 ) lattice artifacts to be smaller and vanish uniformly in 
the 0(a) improved theory. 

Motivated by the apparent linearity of c\ as a function of the lattice spacing within 
our data region, we choose the following interpolation formula for ca as a function of 
5o = 6/0: 


CA(ffo) 


—0.006033 x 


1 + exp 



(4.1) 


which is constrained to reproduce the 1-loop perturbative result of Ref. [12] as g q goes 
to zero. The piece depending exponentially on the bare coupling is inspired by the per¬ 
turbative expression of the lattice spacing in terms of g^, and it is inserted in order to 
capture the linear behavior reflected in Fig. 4. Our non-perturbative results for ca are 
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then compiled in Fig. 5, together with this interpolating function with coefficients 

Po = 9.2056, pi = -13.9847, (4.2) 

which produces a \' 2 /d.o.f. ~ 0.85; they deviate markedly from the 1-loop perturbative 
prediction in the region of our simulations. Eq. (4.1), together with Eq. (4.2), represents 
our final result. This formula can be used at all bare couplings below g q ~ 1.8, together 
with the statistical errors on ca (he., on ca,o as given in Tab. 2), which are ~ 4 — 5% near 
the largest simulated bare couplings and increase to ~ 7 — 8% near the smallest. 

5 Conclusions 

In this work we have determined non-perturbatively ca(s'o), the coefficient required for 
0(a) improvement of axial current matrix elements in lattice QCD with Nf = 3 flavors 
of Wilson quarks, non-perturbative c sw [11] and the tree-level Symanzik-improved gauge 
action. 

The main result is the interpolation formula in Eqs. (4.1) and (4.2), obtained using the 
standard improvement condition, which is imposed together with a variation of boundary 
wavefunctions in the PCAC relation with external states and along a line of constant 
physics in parameter space. This implies that potentially large 0(a) ambiguities in ca are 
avoided and that its remaining intrinsic 0(a) ambiguities disappear smoothly toward the 
continuum limit. Eq. (4.1) for ca(<?o) (with coefficients (4.2)) is valid for bare couplings 
below ~ 1.8 (or, equivalently, for lattice spacings a < 0.09 fm). We have treated the 
topology freezing, encountered in our simulations at the finest lattice spacing (of around 
0.045 fm), by restricting the improvement condition to the trivial topological sector for all 
ensembles considered. 

The gauge held ensembles entering this work were generated with T = 3L/2, in order 
to be re-used for the determination of the axial current renormalization constant Z\, 
see [18] for a preliminary report. After this is completed (the results of which will appear 
in a future publication), axial current matrix elements such as pseudoscalar meson decay 
constants can be calculated precisely from large-volume ensembles of gauge configurations 
at lattice spacings typically employed in the context of phenomenological applications of 
lattice QCD. 
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